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Abstract-  This  paper  aims  to  describe  a  hybrid  technique  that 
combines  the  feasibility  of  our  recently  developed  Multiresolution 
Analysis  of  Signal  Subspace  Invariance  Technique  (MASSIT)  [1] 
with  the  Finite  Element  Method  (FEM)  analytic  model  developed 
in  [2]  to  obtain  accurate  localization  scheme  of  neural  sources  in 
an  extracellular  recording  environment.  The  power  of  the 
proposed  method  stems  from  the  fact  that  robust  array  signal 
processing  approach  is  fused  with  the  FEM  analysis  yielding  the 
closest  scenario  to  practical  experimental  situations.  Results  from 
experimental  signal  and  noise  simulated  composite  are 
summarized  and  the  overall  performance  is  evaluated. 

Keywords  -  Array  processing,  multichannel  neural  probes,  source 
localization. 

I.  Introduction 

The  need  for  accurate  methods  to  determine  neural  source 
location  relative  to  the  recording  micro-device  in  extracellular 
recording  of  neural  activity  in  the  brain  becomes  inevitable. 
Source  localization  helps  to  understand  connectivity  of  small 
populations  of  neural  cells  and  allows  tracking  of  electrode 
motion  without  losing  continuous  recording  of  individual  cells. 
Moreover,  as  cell  localization  schemes  improve,  controlled 
pharmaceutical  delivery  to  sub-circuits  of  neural  populations 
through  special  microprobe  designs  becomes  more  feasible. 

Generally  speaking,  the  amount  of  information  to  be 
revealed  about  neural  cell  location  is  very  limited  because  of 
the  tremendous  challenges  that  face  the  implanted  system. 
These  range  from  the  limited  bandwidth  of  the  telemetry 
system  to  the  limited  computational  power  that  can  be 
expected  from  the  implant  signal  processor. 

The  main  objective  of  this  paper  is  to  describe  a  hybrid 
signal  processing  scheme  that  integrates  the  information 
revealed  by  our  recently  developed  Multiresolution  Analysis  of 
Signal  Subspace  Invariance  Technique  (MASSIT)  for  blind 
source  identification  described  in  [1,3]  and  the  FEM  analysis 
of  neural  tissue-microprobe  interface  environment  summarized 
in  [2],  Each  of  these  techniques  reveals  different  information 
that  can  be  used  to  localize  the  neural  source  and  they  can  be 
considered  to  complement  each  other.  The  integration  of  the 
overall  information  into  a  single  array  signal  processor  should 
be  able  to  achieve  the  required  goal  of  this  work. 

II.  Methodology 

Due  to  the  nature  of  the  neural  signal  with  its  well  time- 
localized  events  or  “spikes”,  the  signal  processing 
methodology  should  rely  on  time  and  frequency  localization 
characteristics.  These  characteristics  are  strongly  met  by  using 
the  Discrete  Wavelet  Transform  (DWT)  as  a  front-end  stage  of 
MASSIT.  The  key  feature  in  MASSIT  is  that  it  fuses  the  well 
developed  techniques  for  multiresolution  analysis  of  the  signal 
with  the  well  developed  theory  of  array  processing  [4]  in  a 
unified  framework.  We  briefly  describe  below  the  main 


approach  of  the  technique  but,  for  lack  of  space,  we’ll  not 
attempt  to  provide  mathematical  details,  which  can  be  found  in 

[1,5]. 

This  section  is  subdivided  in  three  parts.  First,  the  key  idea 
of  MASSIT  is  described.  Next,  source  separation  is 
accomplished  based  on  spatial  filtering  using  signal  subspace 
results  from  MASSIT.  Finally,  the  approach  for  source 
localization  is  linked  to  the  FEM  results  of  [2],  The  overall 
approach  and  preliminary  results  are  summarized  in  the  last 
section. 

A.  Multiresolution  Analysis  of  Signal  Subspace  Invariance 
Technique  (MASSIT) 

The  first  step  of  the  MASSIT  involves  performing  a 
Stationary  Discrete  Wavelet  Packet  Transformation  (SDWPT) 
[3]  to  the  multichannel  observation  data  matrix  denoted 

Y E  91 Mx  ,  where  M  denotes  the  number  of  channels,  and  N 
denotes  the  number  of  snapshots  acquired  by  the  array  in  time 
interval  T.  Assuming  there  are  P  neural  sources  impinging  on 
the  recording  array,  the  classical  model  for  Y  is  expressed  as 
follows: 

Y  =  AS  +  Z  (1) 

where  As  SiAhP  denotes  the  mixing  matrix, 
Se  Si /,v  V  denotes  the  signal  matrix  to  be  estimated,  and 

Ze  Si MlA  denotes  an  iid  additive  noise  component  both 
spatially  and  temporally  correlated.  The  technique  doesn’t 
assume  any  restriction  on  the  cross  correlation  between  the 
signal  and  noise,  due  to  the  strong  background  neural  activity, 
which  represents  a  major  component  of  the  observed  noise 
process.  This  approach  is  the  strongest  and  most  realistic 
among  all  other  array  processing  techniques  [6]  due  to  the 
nature  of  the  neural  signal  environment.  The  obtained 

transformed  multichannel  coefficient  matrix  Y  E  91 U'  ' 
expresses  Y  in  the  wavelet  domain  in  the /h  subband. 

In  a  second  step,  the  technique  performs  advanced  signal 
processing  computations  involving  Singular  Value 
Decomposition  (SVD)  of  the  spatial  covariance 

matrix RJYE  Si‘xhM  of  Y.  [7],  From  a  theoretic  viewpoint,  the 

information  about  the  neural  source’s  spatial  amplitude 
distribution  or  ‘  footprint  ”  on  the  array  side  is  contained  in  this 
covariance  matrix.  The  SVD  allows  one  to  assess  the  P 
principal  sources  impinging  on  the  array  to  obtain  an  estimate 
of  the  original  source  data  matrix  S  at  the  array  interface. 

The  third  step  in  the  algorithm  estimates  the  mixing  matrix 
A  in  each  of  the  wavelet  subbands  from  second  order  statistics 
or  by  estimating  the  signal  and  noise  subspaces  from  the  output 
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of  the  SVD  stage.  The  last  step  consists  of  selecting  certain 
subbands  from  the  full  wavelet  expansion  that  are  guaranteed 
to  span  the  signal  subspace.  The  criterion  by  which  a 
characteristic  wavelet  tree  is  selected  for  each  source  assumes 
that  the  signal  subspace  remains  invariant  across  the 
multiresolution  levels  obtained  in  the  wavelet  domain.  This  is 
guaranteed  to  be  true  only  in  subbands  where  the 
corresponding  wavelet  basis  functions  span  the  signal’s 
functional  space.  The  separation  of  the  P  neural  sources  is 
achieved  using  a  simple  search  in  the  wavelet  decomposition 
tree  to  find  nodes  that  best  describe  the  signal  energy  based  on 
three  components  describing  the  signal’s  energy  distribution  in 
space,  time  and  scale  domains.  These  are: 

1-  The  eigenvalue  distribution  across  tree  nodes, 

2-  The  characteristic  tree  shape. 

3-  The  eigenvector  describing  the  signal  subspace. 

Details  of  the  algorithm  can  be  found  in  [5].  Fig.  1  illustrates  a 
schematic  of  the  system. 

It  is  worthy  to  mention  that,  due  to  the  sparsity  of  the 
neural  signal  on  the  time  base,  a  spike  detection  stage  usually 
precedes  the  MASSIT  stage  so  that  the  N  snapshots  of  the 
array  contain  only  a  single  event.  It  is  likely  that  the  length  N 
time  window  may  contain  more  than  a  single  event  belonging 
to  different  sources.  In  this  latter  case,  the  MASSIT 
automatically  determines  the  number  of  sources  P  within  N 
using  a  statistical  Likelihood  Ratio  Test  (LRT)  and  performs 
the  appropriate  separation  [5,7].  This  case  is  out  of  the  scope  of 
this  paper.  It  will  be  henceforth  assumed  that  the  analysis 
described  in  the  sequel  is  conducted  for  each  detected  event 
separately. 

B.  Spatial  filtering  for  neural  source  separation 

Once  the  mixing  matrix  A  is  obtained  from  MASSIT,  it  is 
feasible  to  obtain  an  estimate  of  the  deterministic  signal  matrix 
S  without  estimating  the  noise  correlation  due  of  the 
compactness  property  of  the  transform  that  allows  thresholding 
small  coefficients,  an  operation  often  referred  to  as  demising 
[6],  Another  way  to  estimate  S  can  use  the  sparsity  of  the 
neural  signal  on  the  time  base  as  stated  earlier  allowing  the 

noise  spatial  covariance  matrix  Rz  G  to  be  estimated 

from  “ spike  free ”  data  instants.  Once  this  is  achieved,  the 


minimum  variance  BLUE  (Best  Linear  Unbiased  Estimator 
[7])  estimate  of  S  can  be  obtained  by  using  the  signal  subspace 
estimate  from  MASSIT.  This  is  be  expressed  as 

S  =  (ATR-jA)-1  AtR‘Y  (2) 

The  term  that  pre-multiplies  the  observation  matrix  Y  is  just 
the  Pseudo  inverse  of  the  mixing  matrix  A  that  is  used  to 
cancel  out  all  the  interference  and  correlated  noise  components 
once  the  principal  sources  in  Y  are  separated.  If  the  noise 
covariance  is  diagonal  (spatially  uncorrelated),  then  this  is 
equivalent  to  applying  appropriate  weights  to  every  channel  to 
maximize  the  SNR.  If  the  noise  is  spatially  colored,  then  the 
inverse  of  the  noise  covariance  acts  as  a  spatial  filter  for  all 
non-principal  components  in  Y. 

C.  Localization  of  neural  sources 

The  MASSIT  relies  on  estimates  of  the  signal  subspace 
from  data  measurements  and  makes  use  of  them  to  identify  the 
best  wavelet  tree  representing  each  source.  The  array  manifold, 
defined  as  all  the  array  responses,  i.e.  steering  vectors  obtained 
as  the  signal  parameters  vary  over  the  entire  parameter  space 
can  be  obtained  by  modeling  the  array  geometry  using  FEM  in 
the  extracellular  tissue  environment.  The  intersection  of  the 
estimated  signal  subspace  for  each  neural  source  and  the 
modeled  array  manifold  will  immediately  enable  the  solution 
to  the  neural  source  localization  problem  in  the  absence  of 
noise. 

The  study  performed  in  [2]  reveals  some  important  results 
about  how  the  array  manifold  can  be  modeled.  By  modeling 
the  tissue  media  through  which  the  signals  pass  from  the 
source  to  the  electrode  array  as  a  homogeneous  resistive 
conductor  with  no  space  charge,  the  FEM  solutions 
demonstrate  the  distortion  of  the  electrical  field  lines  by  the 
presence  of  the  non-conducting  substrate  of  the  array.  The  only 
good  conductors  are  the  recording  sites  at  which  the  signal 
amplitude  is  estimated  by  the  MASSIT  stage  and  they  are  too 
small  to  distort  the  field  further.  Since  it  is  assumed  to  have  a 
poor  interface  to  the  tissue,  the  substrate  interferes  with  the 
flow  of  current  from  a  source  to  a  sink  surrounding  the  field  of 
interest. 


Fig.  1 :  Schematic  of  the  MASSIT:  The  blocks  W ,  W*1  imply  a  SDWPT,  and  inverse  SDWPT  operators  respectively.  There  are  7+1  nodes  in  the 
wavelet  packet  tree  (node  0  being  the  root  level,  i.e.,  the  time  domain  signal,  and  7=  2L*'-2,  where  L  is  the  number  of  resolution  levels,  i.e.,  tree  depth). 
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This  phenomenon  can  be  interpreted  as  an  amplification  of 
the  observed  signal  particularly  when  the  cell  is  directly  over 
the  substrate.  Consequently,  the  signal’s  amplitude  cannot  be 
assumed  to  attenuate  according  to  a  simple  1/r  or  1/r2  roll-off. 
This  implies  that  there  exists  a  point  spread  function  for  the 
source  on  the  substrate  that  depends  on  the  altitude  and  shift 
from  the  substrate  center.  In  otherwords,  the  “footprint”  of  the 
cell  on  the  array  can  be  deconvolved  to  a  point  in  space.  Fig. (2) 
illustrates  the  point  spread  function  hypothesis  of  the  source 
relative  to  its  position  from  the  substrate  center  at  different 
altitudes  in  two  dimensional  space. 
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Fig.  2:  (a)  Approximate  cell  model  in  3  different  positions  X,  Y  and  Z 
above  a  ID  substrate  with  4  sites.  Altitudes  above  substrate  are  0.05, 
0.1  and  0.15  respectively.  The  substrate  width  is  0.2,  and  site  spacing  is 
0.06.  The  origin  is  assumed  at  the  substrate  center.  The  source  radius  is 
0.02  and  the  zero  equipotential  surface  is  at  radius  1  from  the  origin. 
For  clarity,  the  schematic  is  not  to  scale,  (b)  FEM  solution  for  signal 
strength  at  X  (red) ,  Y  (green),  and  Z  (blue)  as  function  of  the  substrate 
shift  from  the  origin. 


The  localization  of  the  neural  source  is  contingent  upon 
estimating  this  point  spread  function  and  coupling  it  with  the 
signal  subspace  estimated  at  the  array. 

The  problem  can  thus  be  envisioned  as  two  fold: 

1-  On  the  probe  computational  side,  the  spatial  covariance 
matrix  can  be  partitioned  to  subtract  out  common  signal 
profiles  of  correlated  sources,  keeping  only  the 
independent  uncorrelated  part  of  each  source  to  be 
localized. 

2-  On  the  brain  side,  the  probe  geometry  is  translated  into 
gain  functions  that  equalize  the  disturbance  of  the  field  at 
the  recording  sites. 


By  carefully  examining  the  profiles  in  Fig.2-b,  less 
skewness  and  higher  peaks  near  the  substrate  center 
characterize  low  altitude  sources.  High  altitude  sources  have 
broader,  skewed  and  have  lower  peaks  near  the  center.  Since 
skewness  is  a  function  of  the  site  position  on  the  substrate, 
each  probe  geometry  has  an  associated  gain  pattern  that  has  to 
be  included  in  the  computations  of  the  signal  subspace  to 
provide  the  localization  information. 

III.  Results 

The  proposed  method  was  implemented  on  a  simulated  data 
set  extracted  from  experimental  data.  A  total  of  P  =  4  neural 
sources  were  detected  across  a  subarray  of  4  channels  on  a  2D 
array  arranged  as  shown  in  Fig.3. 


Fig.  3:  2D-  16  channel  probe  used  for  the  data  acquisition.  Results  for  a 
subarray  of  channels  6,7,  10  and  1 1  were  the  neural  activity  is  highest  are 
shown  in  Fig.  4. 

Template  waveforms  were  obtained  by  averaging  multiple 
realizations  of  these  sources  across  time  [8].  The  spike 
waveforms  are  shown  in  Fig.  4  for  each  channel  of  the 
subarray.  Experimental  noise  from  spike-free  data  was 
extracted  from  the  same  subarray  data. 


Fig.4.:  Template  waveforms  for  4  neural  sources  detected  across  the 
subarray  of  channels  6,  7,  10  and  1 1  of  Fig.  3. 


For  predetermined  signal  to  noise  ratios,  the  extracted  noise 
was  scaled  and  the  template  waveforms  were  added  to  the 
scaled  noise  at  time  intervals  obtained  from  Poisson  processes 
of  known  parameters.  The  simulation  is  thus  the  closest 
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scenario  to  a  real  experimental  data  set  with  the  difference  that 
event  occurrence  times  were  known  before  to  the  noise 
addition.  This  enabled  to  assess  accurate  performance  of  the 
technique  under  different  SNRs.  The  simulated  data  set  was 
processed  with  the  MASSIT  and  the  signal  subspace  estimate 
of  each  source  is  shown  in  Fig. 5  for  SNR  =  12.5  dB. 


Channel  number 


estimates  of  the  array  response  are  derived  from  a  robust  array 
processing  algorithm  that  estimates  the  neural  source’s 
footprint  at  the  substrate.  The  neural  signal  is  estimated  using  a 
linear  unbiased  estimator  that  has  minimum  variance.  The  next 
stage  applies  gain  functions  derived  from  FEM  modeling  of  the 
substrate  geometry  to  correct  for  the  disturbance  of  the  field 
caused  by  the  presence  of  the  nonconductive  substrate.  The 
outcome  is  used  to  estimate  the  neural  cell  location  relative  to 
the  recording  probe. 

Current  work  is  aimed  to  implement  FEM  models  for  each 
probe  geometry  and  verify  the  consistency  of  the  technique 
with  different  array  manifolds.  Meanwhile,  histological 
methods  involving  staining  techniques  are  being  conducted  in 
our  lab  to  accurately  determine  the  probe  location  after  retract 
to  verify  the  accuracy  of  the  technique  in  locating  sources  in 
the  proximity  of  the  probe  trace. 
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Fig.  5:  MASSIT  signal  subspace  estimates  for  SNR  =  12.5  dB 

The  array  manifold  for  the  2D  probe  geometry  of  Fig.  3  can 
be  obtained  using  a  3D  FEM  similar  to  the  linear  array 
illustrated  in  Fig.2-a.  Studies  are  underway  to  integrate  the 
FEM  method  in  3  dimensions  with  histological  findings.  The 
output  is  directly  applied  to  the  signal  solution  obtained  from 
MASSIT  to  adjust  for  the  substrate  presence  in  the 
extracellular  field.  The  parameter  vector  to  be  estimated 
consists  of  the  shift,  depth  and  altitude  from  the  substrate 
center  in  the  x,  y  and  z  directions,  respectively.  The  array 
manifold  for  variable  parameter  vectors  is  calculated.  Each 
parameter  vector  generated  4  profiles  corresponding  to  the  4 
sites  similar  to  each  altitude  in  Fig.2-b.  The  closest 
intersection,  in  a  mean  square  sense,  with  the  estimated  signal 
subspace  from  MASSIT  enables  immediately  to  obtain  the 
parameter  vector  describing  the  localization  information  from 
the  assumed  origin  (substrate  center). 

VI.  Conclusion 

We  have  presented  a  new  methodology  for  solving  the 
localization  problem  of  neural  sources  in  extracellular 
recordings.  The  method  relies  on  a  two-stage  process  where 
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